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We present a new two-dimensional numerical code called Nada designed to solve the full Ein- 
stein equations coupled to the general relativistic hydrodynamics equations. The code is mainly 
intended for studies of self-gravitating accretion disks (or tori) around black holes, although it is 
also suitable for regular spacetimes. Concerning technical aspects the Einstein equations are formu- 
lated and solved in the code using a formulation of the standard 3-1-1 (ADM) system, the so-called 
BSSN approach. A key feature of the code is that derivative terms in the spacetime evolution 
equations are computed using a fourth-order centered finite difference approximation in conjunction 
with the Cartoon method to impose the axisymmetry condition under Cartesian coordinates (the 
choice in Nada), and the puncture/moving puncture approach to carry out black hole evolutions. 
Correspondingly, the general relativistic hydrodynamics equations are written in flux-conservative 
form and solved with high-resolution, shock-capturing schemes. We perform and discuss a number 
of tests to assess the accuracy and expected convergence of the code, namely (single) black hole 
evolutions, shock tubes, and evolutions of both spherical and rotating relativistic stars in equilib- 
rium, the gravitational collapse of a spherical relativistic star leading to the formation of a black 
hole. In addition, paving the way for speciflc applications of the code, we also present results from 
fully general relativistic numerical simulations of a system formed by a black hole surrounded by a 
self-gravitating torus in equilibrium. 

PACS numbers: 04.25.Dm, 04.40.Dg, 04.70.Bw, 95.30.Lz, 97.60.Jd 



I. INTRODUCTION 

Self-gravitating tori (or thick accretion disks) orbit- 
ing black holes (BHs) are common end-products in a 
number of scenarios of relativistic astrophysics. Theo- 
retical evolutionary paths predict that they may form 
after the merger of a binary system formed by a BH 
and a neutron star (NS) or the system formed by two 
NS (see e.g.[lj and references therein). In addition, they 
may also be the result of the gravitational collapse of 
the rotating core of massive stars 0, Q • State-of-the-art 
numerical simulations have started to provide quantita- 
tive estimates of the viability of such systems to form 
[U, i, i, i 0, S i, [M [O, [H, 111, P . Furthermore, such 
a thick disk plus BH system is thought to be the central 
engine for gamma-ray bursts [H, [la, [13, • Therefore, 
understanding the formation and dynamics of such sys- 
tems is a highly relevant enterprise, along with discerning 
its stability properties. Namely, whether they may be 
subject to axisymmetric and non-axisymmetric instabili- 
ties which could also lead to the emission of a significant 
amount of gravitational radiation. 

In particular, the so-called runaway instability, first 
found by Abramowicz, Calvani and Nobili is an ax- 
isymmetric instability that could destroy the torus on 
dynamical timescales. The numerical study of the run- 
away instability in general relativity has so far been in- 
vestigated under different assumptions and approxima- 
tions (see e.g. [1^ [IJ and references therein). Despite 
some progress has been made the existing works are not 
still conclusive on the likelihood of the instability, mainly 
due to the absence of important physics in the modeling. 



The complexity of handling the presence of a spacetime 
singularity in addition to the (magneto-)hydrodynamics 
and the self-gravity of the (possibly magnetized) accre- 
tion torus, turn very challenging the task of carrying out 
full general relativistic simulations of BH surrounded by 
a self-gravitating torus. The code we present in this pa- 
per is built with this mid-term goal in mind. 

Numerical studies of the dynamics of matter around 
BHs are abundant in the literature (see e.g. the refer- 
ences reported in [12]). For the kind of specific work we 
discuss in this paper we note, in particular, that [H, 
already focused on the numerical evolution of matter in 
dynamical axisymmetric BH spacetimes. Recent devel- 
opments in numerical relativity associated with the punc- 
ture method for evolving BHs [1^, [2^ have proved essen- 
tial to perform accurate and long-term stable evolutions 
of spacetimes containing BHs, and the first successful 
simulations of BH/BH binaries have only been possible 
very recently (see ^7] and references therein). Shibata 

[H, 

[l^ used the puncture method to investigate the 
merger of BH/NS binary system in full general relativity, 
and Faber et al. [1^ considered relativistic spherical ac- 
cretion onto a BH. In addition, Baiotti and RezzoUa [1^ 
were also able to handle the accretion of matter onto a 
newly formed BH resulting from the collapse of a rotating 
NS without needing to excise the singularity. 

In order to take advantage of this approach, Shi- 
bata [soj presented a formulation for computing equi- 
librium configurations of a BH and a self-gravitating 
torus in the puncture framework. Previous computa- 
tions of this system by Nishida and Eriguchi [31j, and 
also by Ansorg and Petroff [s^l were not done in the 
puncture framework. In this paper we use the initial 



2 



data in the formalism presented in 30] to carry out 
the first dynamical simulations of a self-gravitating torus 
in equilibrium around a BH as a first step towards a 
systematic investigation of the runaway instability in 
full general relativity. Such initial data are evolved us- 
ing a new two-dimensional, general relativistic hydrody- 
namics code called Nada which is presented and thor- 
oughly tested in this paper. As we explain in detail be- 
low our code solves the 3-1-1 Einstein equations using 
the formulation originally proposed by Nakamura 33 1 
and subsequently modified by Shibata-Nakamura 34 1 
and Baumgarte-Shapiro [35| . which is usually known as 
the BSSN formulation. On the other hand, the gen- 
eral relativistic hydrodynamics equations are written in 
flux-conservative form and solved with a high-resolution 
shock-capturing scheme (HRSC) [2^, |3^ [33|. In addi- 
tion, the Nada code implements the Cartoon method [l^ 
which allows to impose the axisymmetry condition while 
still admitting the use of Cartesian coordinates, the punc- 
ture/moving puncture approach to deal with the BH sin- 
gularity and derivative terms in the spacetime evolution 
equations are computed using a fourth-order centered fi- 
nite difference approximation. 

The code is written in FORTRAN 90, requires approx- 
imately 8 bytes of memory per grid point, and it takes 
about 200 microseconds per grid point per time step us- 
ing a five stages Runge-Kutta (see section III) for the 
time integration on a 2.0 GHz AMD Opteron Dual Core 
270. 



A. Formulation of Einstein equations 

1. BSSN formulation 

We follow the 3-1-1 formulation in which the spacetime 
is foliated into a set of non-intersecting spacelike hyper- 
surfaces. In this approach, the line element is written in 
the following form 



-{a^ - P^f3')dt^ + 2P,dx'dt + j^jdx'dx\ (2.1) 



where a, /3* and 7.^ are the lapse function, the shift three- 
vector, and the three-metric, respectively. The latter is 
defined by 



(2.2) 



where n'^ is a timelike unit-normal vector orthogonal to 
a spacelike hypersurface. 

A reformulation of the ADM system, the BSSN for- 
mulation [33, ^] , has been implemented to solve the 
Einstein equations. Initially, a conformal factor is intro- 
duced, and the conformally related metric is written as 



(2.3) 



such that the determinant of the conformal metric, 7^, is 
unity and (f> — ln(7)/12, where 7 = det(7y). We also de- 
fine the conformally related traceless part of the extrinsic 



curvature K, 



The organization of the paper is as follows: the for- 
mulation of the Einstein equations, including the imple- 
mentation of the puncture approach and gauge condi- 
tions, along with the formulation of the general relativis- 
tic equations is briefly presented in Section II. Section 
III gives a short description of the boundary conditions 
and numerical methods employed for the time evolutions. 
Section IV and V are devoted to present results from tests 
the code has passed, for vacuum and non- vacuum space- 
times, respectively. Section VI discusses numerical simu- 
lations of self-gravitating equilibrium tori, in preparation 
for subsequent work. A summary of our conclusions is 
given in Section VII. We use units in which c = G = 1, 
and in sections V. C to V. F we also use Mq = 1. Greek 
indices run from to 3, Latin indices from 1 to 3, and we 
adopt the standard convention for the summation over 
repeated indices. 



II. BASIC EQUATIONS 



(2.4) 



where K is the trace of the extrinsic curvature. After 
introducing these new variables, the evolution equations 
can be expressed as 



+ laK^+4Tra {E + S) , 



(2.5) 
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(2.6) 
(2.7) 

(2.8) 



-2A'^d,a + 2a ( r]f^A>'' 



+ 1^^0,(3^ + ^j%if3^ + 7'^%/3\ (2.9) 



We give next a brief overview of the formulation for 
the system of Einstein and hydrodynamic equations as 
have been implemented in the code. 



where F* ^j. are the connection coefficients associated 
with 7jj, Cp refers to the Lie derivative along the shift 
vector (see e.g. [1^ for the tensor weights necessary to 
evaluate the Lie derivatives), and Di and i?y are the 
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covariant derivative operator and the three-dimensional 
Ricci tensor associated with the three-metric , respec- 
tively. The object I", known as the conformal connection 
functions, is defined as 

r = j^'^r = -djf^ . (2.10) 

The Ricci tensor Rij that appears in the source term 
of the evolution equation (j2.8p is split into two parts as 
follows 

%=4.-fi?y, (2.11) 

where i?^ is given by 

Rf^ = -2D,Dj(f> - 27yD'=Dfe(/. 

-f 4A0£'j0-47zj£'V^fe0, (2.12) 

where Di is the covariant derivative with respect to the 
conformal metric 7y. The conformal Ricci tensor, Rij, 
is expressed as 

+ 7™" (2f ^;,(,f,),.„ + f f„f fc„„) . (2.13) 

During the evolution we also enforce the constraints 
Ti^Aij) — and det(7y ) = 1 at every time step by using 
the following substitutions 

i,,^i,,--^^7.„ (2.14) 

(2.15) 

The matter source terms, E, Si, and Sij appearing in the 
Einstein equations are projections of the stress-energy 
tensor T'^'^ on the hypersurface with respect to the unit 
normal n'^ 

E = n^n^T'^'', 
5,, = J.^.lj.T'^-, (2.16) 

with s s,jY^. 

In addition to the evolution equations there are three 
constraint equations, the Hamiltonian, the momentum 
and the Gamma constraints, which are only used as di- 
agnostics of the accuracy of the numerical evolutions 
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- —K"^ + 2T:e>'t'E = 0, 


(2.17) 






= 0, 






(2.18) 




= r + djf^ = 0. 


(2.19) 



2. Puncture approach 

Recent breakthroughs in numerical relativity have fi- 
nally made possible accurate and long-term stable 3-1-1 
evolutions of singular spacetimes, including the challeng- 
ing cases of the collision of compact binaries formed by 
either two BHs or a BH and a NS (see [27] and references 
therein) . One of the key ingredients for such success has 
been the so-called "puncture approach" [3§| , in which BH 
initial data are modeled by the Brill-Lindquist topology 
(40l |. where a "throat" at the BH horizon connects two 
asymptotically flat regions. One of the asymptotically 
flat ends is compactifled to a single point known as the 
puncture, leading to a coordinate singularity. 

The puncture approach has the advantage that its nu- 
merical implementation within the BSSN formalism is 
rather simple. The original proposal only considered the 
fixed puncture approach (4ll . |43 . |43| , where the confor- 
mal factor is split into a regular part and a singular part. 
Although this method does not lead to long-term sta- 
ble evolutions of spacetimes containing BHs, it allows for 
a number of code tests. However, two different groups 
[25l.[26j developed recently the moving puncture approach 
in which no singular term of the conformal factor is fac- 
tored out, and the punctures are allowed to move through 
the grid. The difference between these two methods is 
how the conformal factor is evolved. In the so-called 
(/>- method [1^, the original BSSN variable 4> is evolved 
through the usual BSSN evolution Eq. (|23)) . On the 
other hand, the so-called x-method 25] introduces a new 
conformal factor deflned as x = e"**"^, and the following 
evolution equation, that replaces Eq. (|2.5p . 

{dt~Cp)x^\x{o^K-d,P^). (2.20) 

This moving puncture approach, together with the so- 
called puncture gauge (see Sec. II A. 3), has led to ma- 
jor success in simulations of binary BHs (H, [i^, 113, 
S, m, HO]. Motivated by this, we have implemented 
both methods in our 2D code, the (/)-method and the x~ 
method, for the moving puncture approach to investigate 
the dynamics of matter around BHs. 

3. Gauge choices 

In addition to the BSSN spacetime variables 
{•^ij,Aij,K,(l) or Xj^*); there are two more variables 
left undetermined, the lapse, a and the shift vector, (3^. 
The code can handle arbitrary gauge conditions, however 
long-term BH evolutions in the moving puncture frame- 
work have been successful with some combination of the 
so called "1-l-log" condition ^ for the lapse, and the 
"Gamma- freezing" condition for the shift vector [43*] . 

The form of this slicing condition most commonly used 
for moving puncture evolutions contains the advective 
term [3^diOi and is expressed as 

dta - P'dia = -2aK. (2.21) 
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Several authors (see, e.g. [2^, 153, l53| ) have used the so- 
called "non-advective l+log", by dropping the advective 
term in Eq. (j2.2ip . for single BH evolutions. In this case, 
the slicing condition takes the form 



dta ~ —2aK. 



(2.22) 



It was shown by [52| , that the evolution of a single punc- 
ture using the condition given by Eq. (|2.22p settles down 
into a time-independent and maximally sliced solution. 

For the simulations of a BH surrounded by a self- 
gravitating torus presented in this paper, we choose the 
slicing given by Eq. (|2.22p as the BH does not move 
through the grid, and therefore the advective term (3'^dia 
does not play an important role. 

For the shift vector, we choose the Gamma-freezing 
condition, written as 



(2.23) 



iiW 



Ui = hui, 



pb<P p 

e=—T^,n^^n-' = hW-—, 
p* pW 



(2.27) 
(2.28) 

(2.29) 



W = au\ (2.30) 

where W and h are the Lorentz factor and the specific 
fluid enthalpy respectively, and P is the pressure. 

By defining the vector of unknowns, U, and fluxes 
and F^ along the x and z directions as 



(2.31) 



+ PaV7(«"-h /?")], (2.32) 



dtB' = dtV' - rjB\ (2.24) 

where 77 is a constant that acts as a damping term, in- 
troduced both to prevent long term drift of the metric 
functions and to prevent oscillations of the shift vector 
[4^. This parameter also has an effect on the coordi- 
nates of the final spacelike hypersurface. In Ref. [s^ it 
was shown that increasing the value of 77 increases the co- 
ordinate size of the BH, which allows for better numerical 
resolution across the BH. On the other hand, larger val- 
ues of the damping parameter introduce a higher drift in 
the location of the horizon in time, as well as in the de- 
formation of the metric during the evolution, caused by 
the larg er values reached by the connection functions F' 
(28I Is^f. Bearing this in mind, we use 77 — 0.3/Af, where 
M is the ADM mass of the system, for the evolutions of 
spacetimes containing a BH presented in this paper. 

B. Formulation of the hydrodynamics equations 

The general relativistic hydrodynamics equations, ex- 
pressed through the conservation equations for the stress- 
energy tensor T^'' and the continuity equation are 

V^Tf^" - , (p?i^) = 0, (2.25) 

where p is the rest-mass density, is the fluid four- 
velocity and V is the covariant derivative with respect to 
the spacetime metric. Following [37| the general relativis- 
tic hydrodynamic equations are written in a conservative 
form in cylindrical coordinates. Since the Einstein equa- 
tions are solved only in the y = plane with Cartesian 
coordinates, the hydrodynamic equations are rewritten 
in the Cartesian coordinates for y = 0. The following 
definitions for the hydrodynamical variables are used 

p^ = pWe^"^, (2.26) 



F^ = J^v\ Jyv\ J^v'- + Pay^, 

E,v' + Pa^iv' + P')], (2.33) 

with Ji = p.^,xli and = p^e, the set of hydrodynamic 
equations (|2.25p can be written in conservative form as 

dtU + a^^F^ + a^F^ = S, (2.34) 

where S is the vector of sources. We refer to [37*1 for 
further details on these equations, in particular regarding 
the form of the source terms. 

To close the system of equations, we choose two pos- 
sible equations of state, the so-called F-law equation of 
state (ideal fluid) given by 

P=(r-l)pe, (2.35) 

where e is the specific internal energy, and a polytropic 
equation of state 

P = np^. (2.36) 

Here k is the polytropic constant, F = 1 -I- and 
N is the polytropic index. In those simulations where 
the system evolves adiabatically, such that no shocks are 
present, we use the polytropic equation of state during 
the evolution. 

After each time iteration the conserved variables 
(i.e. p*, Jj;, Jy, Jz, i?*) are updated and the primitive hy- 
drodynamical variables (i.e. p, w^, w^, w^, e) have to be re- 
covered at the corresponding step. This is done by solv- 
ing the following equation for the Lorentz factor, W , de- 
rived from the normalization of the 4- velocity of the fluid 

I^^^l + 7-^.,^, = l + 7^^u.u,(^ + ^) . 

(2.37) 

Once solved for W, the other variables, p,v'^,P,e and h 
are computed from Eqs. (I2.26p - (l2.30p and the equation 
of state. 
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III. NUMERICS 

There are two schemes implemented for the update of 
the numerical solution with time: an iterative Crank- 
Nicholson (ICN) scheme taking two corrector steps [ssj . 
and an optimal strong stability-preserving (SSP) Runge- 
Kutta of fourth-order algorithm with 5 stages (RK4) . 
The RK4 scheme is used for simulations of spacetimes 
containing a singularity, where a high-order of accuracy 
is an important issue. 

We use second-order slope limiter reconstruction 
schemes (both minmod and MC are implemented in the 
code) to obtain the left and right states of the primi- 
tive variables (i.e. p,v^ ,v'^ ,€) at each cell interface, 
and these reconstructed variables are then used to com- 
pute the left and right states of the evolved quantities 
(p*, Jx, Jy, Jz, E^). Next, we use HLLE or Roe approx- 
imate solvers to compute the numerical fluxes in the x 
and z directions. Details on such high-resolution shock- 
capturing schemes are available elsewhere (see e.g. p2| 
and references therein). 

Derivative terms in the spacetime evolution equations 
are represented by second or fourth order centered fi- 
nite difference approximation in a uniform Cartesian 
grid except for the advcction terms (terms formally like 
(3^diu), for which an upwind scheme is used. 

A. Discretization of Axisymmetric Systems: 
"Cartoon" Method 

As pointed out by [11] any given system of equations 
in 3D possessing a rotation symmetry with respect to 
the z-axis can be finite differenced and solved in the 
X — z {y ~ 0) plane alone because of the symmetry con- 
dition. When the system of equations is expressed in 
Cartesian coordinates, partial derivative terms with re- 
spect to the y-coordinatc will appear and these need to 
be computed using the computational domain. For in- 
stance, in a second-order centered difference approxima- 
tion, the values of the derivative in the y-direction of a 
quantity /(x, 0, z), are computed using the values of / at 
the nearest two grid points i.e. f(x, — Ay, z), f(x, Ay, z). 
Then, because of axisymmetry, the y-derivatives can be 
determined in the x — z plane from information contained 
in this same plane. The Cartoon method obtains the 
boundary conditions at y = ±Ay that are necessary to 
evaluate the derivatives in the y-direction by means of a 
rotation about the symmetry axis of the different tensor 
quantities. On the other hand, in a fourth-order centered 
difference approximation, the values of the derivative in 
the y-direction of a quantity f{x, 0, z), are computed us- 
ing the values of / at the grid points /(x, ±2Ay, z) and 
/(a;,±Ay,z). 

The values of the variables at the positions 
(\/ x"^ + Ay^, 0, z) and (■\/ x"^ + 4Ay^, 0, z) are interpo- 
lated from the neighboring grid points using Lagrange 
polynomial interpolation (see e.g. [58| ) with interpolat- 



ing polynomials of degree 2 and 4, depending on the order 
of the finite difference approximation. 



B. Boundary Conditions 

The computational domain is defined as < a; < L and 
< z < L, where L refers to the location of the outer 
boundaries. For simulations with the puncture method 
we used a staggered Cartesian grid to avoid that the lo- 
cation of the puncture at the origin coincides with a grid 
point. A number of different boundary conditions are im- 
plemented in the code. These are imposed for the space- 
time variables or the hydrodynamical primitive variables 
at the inner and outer boundaries as follows: for both 
the spacetime and hydrodynamical variables, 7r-rotation 
symmetry is imposed around the z-axis, and equatorial 
plane symmetry with respect to the z — Q plane. At 
the outer boundaries we impose radiative boundary con- 
ditions [i^. Note that we do not apply this boundary 
condition to the conformal connection functions, F*, for 
which we used static boundary conditions. 



C. Atmosphere treatment 

An important ingredient in numerical simulations 
based on finite difference schemes to solve the hydro- 
dynamic equations is the treatment of vacuum regions. 
The standard approach is to add an atmosphere of very 
low density filling these regions [s^. We follow this ap- 
proach and treat the atmosphere as a perfect fluid with 
a rest-mass density several orders of magnitude smaller 
than that of the bulk matter. The hydrodynamic equa- 
tions are solved in the atmosphere region in the same 
way that is done for the region of the bulk matter. If 
the conservative variables or fall below some min- 
imum value then the values of the conserved quantities 
are set to the atmosphere value. Similarly in the rou- 
tine that converts primitive variables from conservatives, 
if the rest-mass density p or specific internal energy e 
fall below the value set for the atmosphere, this point is 
reset to have the atmosphere value of the primitive vari- 
ables. In particular for simulations of relativistic stars, 
and systems composed of a BH plus a self-gravitating 
torus system, the atmosphere density is usually taken to 
be about 6-8 orders of magnitude smaller than the initial 
maximum rest-mass density. 



D. Diagnostics 

To check the accuracy of the numerical simulations we 
monitor the violation of the Hamiltonian constraint, and 
the conservation of the total rest-mass M*, the ADM 
mass M, and the angular momentum J. We compute 
these various quantities in the y = plane as 
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In axisymmetry the apparent horizon (AH hereafter) 
equation becomes a nonhnear ordinary differential equa- 
tion (ODE) for the AH shape function, h = h{e) 
We have implemented an AH finder that solves this ODE 
by a shooting method using that dgh{9 = 0) = and 
dgh{9 = tt/2) = as boundary conditions. We define the 
mass of the AH as 

M.H = (3.4) 
where A is the area of the AH. 



IV. VACUUM TESTS: SINGLE BH 
EVOLUTIONS 

The Schwarzschild metric in isotropic coordinates is 
used as initial data to test the ability of the code to evolve 
BH spacetimes within the fixed and moving puncture ap- 
proaches. These initial data are such that the 3-metric 
is written in Cartesian coordinates as 

da^ ^ij^idx^ + dy^ + dz^), (4.1) 

where the conformal factor is i/" = (1 + M/2r), M being 
the mass of the BH. Here r is the isotropic radius, = 
x"^ + y'^ + . Thus, the spatial metric takes the form 
7ij = '0^7ij , where the conformal metric is the fiat metric. 
Initially the extrinsic curvature is Kij — 0. 

A. Schwarzschild BH with geodesic slicing 

Following the fixed puncture approach, we first evolve 
these initial data with geodesic slicing, that is setting 
a = 1 and = 0. Although, in these coordinates, the 
numerical evolution is known to last very short time and 
is expected to crash at t = ttM when the spacelike hy- 
persurface reaches the physical singularity, there is an 
analytic solution for the evolution of the spacetime which 
can be used to compare to the numerical solution. 

The top- left panel in Fig. [T] shows the evolution of the 
radial metric component 7^^, that is obtained from the 



Cartesian metric functions, along the diagonal for differ- 
ent resolutions together with the analytic solution at sev- 
eral time steps, performed with second-order finite differ- 
ences in space. The convergence of the code can already 
be seen in this panel, but this is better appreciated in the 
bottom-left panel, in which the Hamiltonian constraint 
is plotted at t = 2M for the three different resolutions. 
Violations of the Hamiltonian constraint are a measure of 
the numerical error, and this figure shows that this error 
scales at the expected rate for second-order convergence. 
Results obtained with a fourth-order finite differencing 
are shown in the upper and lower right panels, where 
the same quantities are plotted using this higher order 
differencing. Clearly, using fourth-order finite differenc- 
ing increases the accuracy of the numerical evolutions, 
although the convergence rate for the highest resolution 
is not exactly fourth-order in the whole computational 
domain. We note that in addition to the numerical error 
due to the finite difference approximations for the spatial 
derivatives, there is another source of numerical error due 
to the interpolation needed with the Cartoon method. 

B. Schwarzschild BH with 1 -f log slicing 

Next, we use a zero shift and a modified 1 -f log slicing 
condition that permits the increase of the life of the sim- 
ulation up to about 30-40Af. The slicing condition used 
for this simulation is [1^, 

dtOL = -2aiP'^K. (4.2) 

Figure [2] displays the time evolution of one of the com- 
ponents of the conformal three metric and the lapse func- 
tion, and it proves how the new gauge condition enables 
for longer simulations. For this simulation we consider a 
grid resolution of Aa; = 0.05A/ with A^^j, x iV^ = 200 x 200 
grid points. As expected, the metric function grows 
due to the grid stretching related to the collapse of the 
lapse. In order to increase further the length of the sim- 
ulations of a single BH spacetime, we have implemented 
the moving puncture method, and its test is shown next. 

C. Schwarzschild BH and the moving puncture 
approach 

It has been pointed out [H, [H, [11], that using the 
"moving puncture method" the numerical slices of a 
Schwarzchild BH spacetime reach a stationary state after 
about AQM of evolution. 52] computed numerically this 
time-independent puncture data for the Schwarzschild 
spacetime and showed that, as the evolution proceeds, 
the initial slice settles down to a maximal slice with 
K — Q. This time- independent data has recently been 
obtained analytically by Baumgarte and Naculich [g^] . 

With the aim of testing the ability of our code to per- 
form long-term, stable and accurate, evolutions of a sin- 
gle BH, we first check the evolution of the BH punc- 
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FIG. 1; Top-left panel: evolution of the conformal metric component 7^,- along the diagonal of the grid. Different lines refer 
to the numerical solution with resolutions of 100, 200, and 400 grid points in each direction. The analytical solution is also 
shown. Bottom-left panel: Hamiltonian constraints along the diagonal at time t = 2M for three different resolutions showing 
the expected second-order convergence. The right panels show the same quantities and the expected convergence rate when 
using fourth-order finite differencing. 
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FIG. 2: Time evolution of the conformal metric component 
and the lapse a along the a;-axis using 1 + log slicing. 

ture analytical data given by (6^ . with the gauge condi- 
tions given by Eqs. (|2.22p - (|2.24p . and with two different 
methods for the evolution of the conformal factor, the 0- 



method, which uses Eq. (|2.5p . and the x-method which, 
instead uses the evolution Eq. (|2.20p . For the damping 
parameter in the evolution equation for the shift, we set 
ri = 0.3/M. 

In the upper panel of Fig. [31 we plot the absolute de- 
viations for the conformal factor between our numeri- 
cal solution and the exact solution sA t = 9M obtained 
with the 0-method, for three runs with different reso- 
lutions (Ax = 0.12M,0.06M,0.03M, with N^x = 
150 X 150, 300 X 300, 600 x 600). Deviations for the con- 
formal factor are shown in the inner region of the grid, 
which at < = 9M is not affected by the outer boundary 
conditions. The inset, shows these deviations from the 
exact solution, but rescaled assuming fourth-order con- 
vergence. We note approximate fourth-order convergence 
for the two lower resolutions runs (Aa; = 0.12Af, 0.06M), 
but, although still converging, the convergence rate de- 
creases from fourth-order for the highest resolution run 
[Ax — 0.03M) particularly in the region outside the AH. 
For this resolution the error between the numerical and 
the analytical solution in the region outside the AH is 
of the order of 10~^-10~^, and we find that increasing 
further the resolution does not reduce significantly this 
error. 

In the bottom panel of Fig. [3] we plot again the devi- 
ations for the conformal factor, however these are now 
computed using the ^-method as well as the 0-method, 
and obtained with the same resolution, Ax = 0.06Af, at 
t = 9Af. Clearly, the deviations from the exact solution 
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FIG. 3: In the top panel we consider as initial data the ana- 
lytical solution proposed by Baumgarte and Naculich 62;3 and 
show the deviations for the conformal factor from the exact 
solution for different resolutions using the (^-method. In the 
inset, we show deviations rescaled to show the convergence. 
In the bottom panel, we show the deviations of the confor- 
mal factor using the same resolution, Aa; = 0.06M, with the 
0-method as well as the x-method. 



obtained with the ^-method are smaller than with the 
(/)-method. In fact, these results are to be expected, as 
(j) diverges like log r near the puncture, while x behaves 
like near the puncture. 

We next perform a comparison between the time- 
independent analytical data [g^I and the late time nu- 
merical solution of a Schwarzschild spacetime time, in 
isotropic coordinates, choosing initially a = ip~'^. For 
this comparison we use the same gauge conditions as 
above, grid spacing Aa; = Az — O.IM and a grid with 
Nx X Nz — 300 X 300 points to cover the computational 
domain. In Fig. |4l we plot with a solid line the time- 
independent analytical solution for the conformal factor, 
while solid circles indicate the numerical solution for the 
conformal factor of the Schwarzschild spacetime initial 
data at i = 200M, computed using the X" method and 
empty circles correspond to the same quantity but us- 
ing the 0-method. Both simulations are compared at a 
sufficiently late time to allow us to evaluate the effect 
of the outer boundary conditions on the late time evo- 
lution. Clearly, the numerical solution agrees well with 
the analytical solution, both for the x and ((> methods. 
However, again the x-method gives better results in the 
region near the puncture. In addition, we observe that 
the outer boundary conditions do not affect the long-term 
stability and accuracy of these simulations. 



HYDRODYNAMIC TESTS 



Relativistic shock tube test 
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B&N-solution, Ax=0.1M 
t = 200M, ;(:-method 
t = 200M, 0-method 




FIG. 4: We show, with a solid line, the conformal factor given 
by Baumgarte and Naculich 63], and with solid circles the 
numerical solution for the conformal factor at t — 200M of 
Schwarzschild spacetime in isotropic coordinates initial data 
with the x-method and with empty circles the same quantity 
but computed using the (/!>-method. 



As a first test of the solution of the relativistic hydrody- 
namic equations we perform a standard one dimensional 
shock-tube problem (a Riemann problem) in flat space- 
time. This is a common test to assess hydrodynamical 
codes (see e.g. [M])- In this test, two uniform and dif- 
ferent fluid states are initially separated by an interface 
which is then instantaneously removed. The initial data 
for the left state is Pl = 13.33, pl = 10.00, = 0.00 
and for the right state Pr = 0.66 x 10~^, pR — 1.00, vr = 
0.00, with the initial interface located at z — 0.5. The 
fluid is assumed to be an ideal fluid with F = 5/3. 

The results of the numerical evolution are shown in 
the left panels of Fig. [5l where we plot the profiles at 
time t — 0.3 of the pressure, density and z-component 
of the fluid velocity for a shock-tube test along the z- 
direction. The solid line represents the exact solution 
of the shock-tube problem computed using the public 
domain code RIEMANN developed by J.M. Marti and E. 
Miiller [g^]. The numerical solution is represented by 
crosses. We use 400 grid points in the z-direction and a 
grid spacing of Az = 1/400. The particular combination 
of schemes used for this hydrodynamical test comprise 
the Roe solver, MC cell-reconstruction, and ICN for the 
time update. 
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FIG. 5: Left panel: comparison of the numerical profiles of density, pressure and velocity with the analytical solutions for the 
solution of a one dimensional shock tube problem at time t — 0.3. Right panel: wall shock problem with v — 0.9 at time 
t = 1.6. The exact solution is represented by the solid line and numerical solution by crosses. 



B. Relativistic planar shock reflection test 

In our second test in Minkowski spacetime we carry out 
a relativistic planar shock reflection problem along the z- 
direction with the following initial conditions P = 10~^, 
p = 1.0 in the entire domain, while v = 0.9 for z < 0.5, 
and V — —0.9 for z > 0.5. We use an ideal fluid equation 
of state with T = 4/3. We show results of the evolution 
in the right panels of Fig.[5l where we plot the profiles of 
the solution at time t = 1.6 for the pressure, density and 
z-component of the fluid velocity. Again, the solid lines 
refer to the analytic solution while the numerical solution 
is represented by crosses. The schemes used for this sim- 
ulation are the HLLE solver, minmod cell-reconstruction, 
and ICN for the time update. 



C. Relativistic spherical shock reflection test 

The initial configuration for this test problem consists 
of a medium with uniform density (p — 1) and pressure 
(P = 2.29 X 10-5(r - 1)) where T = 4/3, with constant 
spherical inflow velocity Vin — —0.9 and Lorentz factor 
Win- Initially, at t = 0, the gas collides at the center of 
symmetry which forms a strong shock wave that propa- 
gates upstream. The analytic solution for this test has 
the following form [g^ 



P{r,t)^{ } , \A (5.1) 

with the compression ratio given by 

^=^ + ^(m„-l), (5.2) 
and the shock velocity by 

Vs = ^WM. (5.3) 

This test problem has been proposed by [i^ to test 
the ability of a three-dimensional special relativistic hy- 
drodynamics code in Cartesian coordinates to keep the 
spherical symmetry of the solution. For our axisymmet- 
ric hydrodynamics code in Cartesian coordinates, this is 
a two-dimensional test problem, which allow us to eval- 
uate the solution of the hydrodynamic equations both in 
the X and z directions. We use x = 400 x 400 grid 
points with Ax = Az = 1/400 to cover the computa- 
tional domain. The schemes used for this simulation are 
the HLLE solver, minmod cell-reconstruction, and ICN 
for the time update. 

We show results of the evolution in Fig. [SI where we 
plot the profiles of the solution at time i = 3 for the 
density, pressure and fluid velocity along the diagonal. 
The solid line represents the analytic solution while the 
numerical solution is represented by crosses. 
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FIG. 6: Comparison of the numerical profiles of density, pres- 
sure and velocity with the analytical solutions for the solution 
of a relativistic spherical shock reflection test at time t — 3 



To quantify the error we computed the relative global 
error between the numerical and the analytical solu- 
tion defined by eroi = Cabs/Ei^fc |w(a;i, z^, i„)|AxAz], 
where the absolute global error is eabs — J2i k l^ffe ~ 
w{xi, Zk,tn)\^xAz, w^f. is the numerical solution and 
w{xi, Z]~,tn) is the analytical solution. We find that the 
relative global errors for the density, pressure and velocity 
are 2.1%, 1.1% and 0.6% respectively. As for the previ- 
ous shock-tube test and planar shock reflection test, the 
agreement achieved between the analytic and numerical 
solution is remarkable. 



D. Spherical relativistic stars in a fixed spacetime 

For the next test we use the solution of the Tolman- 
Oppenheimer-Volkoff equations to assess the capability 
of the code to perform long-term stable numerical sim- 
ulation of a NS in equilibrium. In order to check the 
hydrodynamical evolution independently from the space- 
time evolution we follow the common approach of keeping 
the spacetime fixed during the numerical evolution. This 
is known as the Cowling approximation [TOj, in analogy 
with the corresponding approximation in perturbative 
studies of stellar oscillations in which the metric com- 
ponents (or the gravitational potential in the Newtonian 
framework) are kept constant. 

It has been shown [5^ that the truncation errors of the 
finite difference representation of the PDEs are enough 
to excite small periodic radial oscillations which manifest 
themselves as periodic variations of the hydrodynamical 
(and spacetime) quantities. The power spectrum of the 



evolution of the central density provides the frequencies 
of the fundamental mode of oscillation and of its over- 
tones, which can be compared with the corresponding 
frequencies computed by perturbative techniques. Al- 
though the excited oscillations are purely numerical in 
origin (i.e. their amplitude converges to zero as the reso- 
lution increases) they still represent the oscillation modes 
of the relativistic star and their frequencies should agree 
with the eigenfrequencies computed by linear perturba- 
tion analysis. For the purposes of further testing the 
code and compare with independent results we focus on 
an initial TOV model that has been extensively investi- 
gated numerically by [H, [6^ . This model is a relativis- 
tic star with = 1, polytropic constant k = 100 and 
central rest-mass density pc = 1.28 x 10"'^ so that its 
mass is M = 1.4, its baryon rest-mass = 1.5 and 
its radius R = 9.59. We evolve these initial data with 
our non-linear code and compute the power spectrum of 
the evolution of the central rest-mass density. We note 
that in the simulations of spherical stars in equilibrium 
presented below, we use the polytropic equation of state 
during the evolution. 

In the upper-left panel of Fig. [7] we plot the time evo- 
lution of the central rest-mass density for a simulation 
with Ax = 0.15. We observed that the truncation er- 
rors at this resolution are enough to excite small peri- 
odic radial oscillations, visible in this plot as periodic 
variations of the central density. We see that the damp- 
ing of the periodic oscillations of the central rest-mass 
density is very small during the whole evolution, which 
highlights the low numerical viscosity of the schemes im- 
plemented. By computing the Fourier transform of the 
time evolution of the central rest-mass density we obtain 
the power spectrum, which is shown with a solid line in 
the lower-left panel of the figure, while the dashed ver- 
tical lines indicate the fundamental frequency and the 
first five overtones computed by [gs] with a non-linear 
hydrodynamics code using spherical polar coordinates. 
The agreement found for the fundamental frequency and 
overtones is very good, with the relative error between 
the fundamental frequencies being less than 1%. 



E. Spherical relativistic stars in a dynamical 
spacetime 

For our first test of the coupling of the Einstein equa- 
tions and the general relativistic hydrodynamic equations 
we again use the TOV solution. In analogy to what hap- 
pens in the Cowling approximation, truncation errors ex- 
cite small oscillations in the star. Now, however, the 
truncation errors come not only from the hydrodynamic 
part of the code but also from the spacetime part solving 
the full set of Einstein equations. 

We computed the eigenfrequencies of the coupled evo- 
lution of the same TOV star with k = 100, A^ = 1 and 
central rest-mass density pc = 1.28 x 10"'^ which has been 
discussed above under the assumption of a fixed back- 
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FIG. 7: Top panels show the time evolution of the normalized central density, in the Cowling approximation (left panel), and 
in a dynamical spacetime evolution (right panel). Power spectrum of the evolution of the central rest-mass density for an 
M = 1.4, K = 100, A'^ = 1 polytrope in the Cowling approximation (bottom-left panel). F represents the frequency of the 
fundamental mode and H1-H6 are the first six overtones computed by p3^ . Correspondingly, the bottom-right panel shows the 
same quantities in a dynamical spacetime evolution of the same TOV star. The frequencies computed by ^] are displayed 
with dashed vertical lines. 



ground spacetime. Again, the radial oscillations excited 
by truncation errors manifest in the time evolution of the 
central rest-mass density. The coupling to the spacetime 
increases the amplitude of the oscillations, and also shifts 
the frequencies of the modes towards lower frequencies. 
This is shown in the lower-right panel of Fig. [7] which dis- 
plays the power spectrum of the evolution of the central 
rest-mass density (solid line), and the eigenfrequencies 
computed by [501 (dashed lines). We note that the loca- 
tions of the frequency peaks for the fundamental mode 
and the two overtones is in very good agreement, with 
the relative error in the fundamental frequencies being 
less than 1%. 



F. Gravitational collapse of marginally stable 
spherical relativistic stars 

We next test the capability of the code to follow BH 
formation with the gravitational collapse to a BH of a 
marginally stable spherical relativistic star. For this test, 
we consider a. k — 100, N = 1 polytropic star with central 
rest-mass density pc = 3.15 x lO"'^, its mass is M — 1.64 
and its baryon rest-mass Af* = 1.79. In order to induce 
the collapse of the star we initially increase the rest-mass 
density by 0.5%. 

We present numerical results for the simulation of 
the gravitational collapse of a marginally stable spher- 
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FIG. 8: Time evolution of the normalized central density, 
the central lapse function, and mass of the AH in units of 
the ADM mass of the system for the collapse of a marginally 
stable spherical star to a BH. Dotted lines represent results 
obtained with the lowest resolution, dashed lines with the 
medium resolution while solid lines with the highest resolu- 
tion. 
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ical relativistic star performed with Ax = Az = 
O.IM, 0.075M,0.05M to show the convergence of the 
code. With these resolutions the star is initially covered 
by approximately 60, 80 and 120 points. The hydrody- 
namics equations are solved using the Roe solver and MC 
reconstruction, while the Einstein equations are solved 
with fourth-order finite differencing and the gauge con- 
ditions given in Eqs. with 77 = 0.5. Time 
integration is done with RK4. We plot in Fig.[8]the time 
evolution of the normalized central density (top panel), 
of the lapse function at the center (middle panel), and 
of the mass of the AH (Eq. [33) in units of the ADM 
mass of the system (bottom panel). Dotted lines repre- 
sent results obtained with the lowest resolution, dashed 
lines with the medium resolution while solid lines with 
the highest resolution. Overall, as the collapse proceeds 
the star increases its compactness, and is reflected in the 
increase of the central density as shown in the top panel. 
The middle panel shows the characteristic collapse of the 
lapse function at the center of the star indicating the for- 
mation of a BH. The most unambiguous signature of the 
formation of a BH during the simulation is the forma- 
tion of an AH. Once an AH is found by the AH finder, 
we monitor the evolution of the AH area, and also of 
its mass which is plotted, in the bottom panel of Fig. [H 
This panel shows that the formation of the BH delays 
with decreasing resolution, since the increase of the cen- 
tral density slows down because numerical dissipation is 
larger for the lower resolution runs. An AH is first found 
at approximately eX t = 180Af {t = 0.88 ms) for the sim- 
ulation with Aa; = Az = 0.05M. This figure also shows 
that the mass of the AH relaxes to the ADM mass of the 
system. The difference in the ADM mass and the mass of 
the AH when we stop the evolution at t = 300Af , about 
120M after the AH is first found, is less than 1%. 

This test shows that the code is able to follow BH for- 
mation and its subsequent evolution for many timescales. 
We note that, unlike the results discussed in [2^, we do 
not need to add any numerical dissipation to the evo- 
lution equations for the spacetime variables and gauge 
quantities to perform this simulation. Instead, we rely 
only on the gauge choice used here to follow the forma- 
tion and evolution of the BH formed as a result of the 
gravitational collapse. 



G. Rapidly rotating relativistic stars 

The evolution of stable rapidly rotating relativistic 
stars is a more demanding test than all previous ones, as 
this involves testing parts of the code that are now used 
since there is a non-zero y-component of the shift vector. 
The initial data used for this test are the numerical solu- 
tion of a stationary and axisymmetric equilibrium model 
of a rapidly and uniformly rotating relativistic star with 
angular velocity which has been computed using the 
Whisky code [2^,|71|. The initial data for this equilibrium 
rotating star are initially computed in spherical polar co- 
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FIG. 9: Profile of the y-component of the velocity of a rapidly 
rotating relativistic star along the a;-axis, at the initial time 
(solid line) and at time t = 2.46 ms, which corresponds to 2 
rotational periods. 



ordinates and then transformed to Cartesian coordinates 
using the standard coordinate transformation. 

Here, we consider a uniformly rotating polytropic star 
with r = 2, pc = 1-28 X 10-3, and rotating at 92% of 
the allowed mass-shedding limit for a star with the same 
central rest-mass density. The ratio of the polar to equa- 
torial coordinate radii for this model is 0.7, its mass is 
M — 1.57 and its baryon rest-mass — 1.69. For this 
simulation we use the "1-l-log" condition for the lapse, 
and the "Gamma-freezing" condition for the shift vec- 
tor, with rj — 3.0. We use the so-called F-law equation of 
state during the evolution. 

For this test, the outer boundary of the grid is placed 
about three times the equatorial radius. Wc check that 
during the evolution, the profiles of the different vari- 
ables are kept very close to their initial value for several 
rotational periods, which is shown in Fig. (S] Here we 
display the profile of the y-component of the velocity, , 
along the x-axis at the initial time, t — (solid line) 
and at t = 2.46 ms (dashed line), with the latter corre- 
sponding to two rotational periods. As has been shown 
by [H, , the small difference originating at the surface 
of the star after two rotational periods is expected since 
we used a second order reconstruction method with the 
MC limiter. Nevertheless, it shows that deviations of the 
numerical solution from the initial profile are very small. 
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FIG. 10: Isocontours of the logarithm of the rest-mass density of the torus. The left panel shows the configuration at the initial 
time and the right panel the corresponding distribution after 5 dynamical timescales (approximately lOOOM). The equilibrium 
solution is preserved to high accuracy. 



VI. SIMULATIONS OF SELF-GRAVITATING 
TORUS IN EQUILIBRIUM AROUND A BH 

A. Initial data for BH and self-gravitating torus 
system 

The initial data for the numerical simulations of the 
system formed by a BH and a self-gravitating torus in 
axisymmetric equilibrium have been recently computed 
by [30|. The basic equations for the metric functions 
are derived assuming the 3-1-1 formalism, and the line 
element is written in the quasi-isotropic form as 

ds^ = -a^dt^ +TP^[e^''{dr^ +r^d0^)+r^ sin^ 0{l3dt+dip)^]. 

(6.1) 

Where e^^ denotes the conformal metric for the rr and 69 
parts. The equations for the metric functions are solved 
in the BH puncture framework, and it is assumed that 
the puncture is located at the origin. 

The equilibrium configurations for the matter are ob- 
tained by assuming a perfect fluid stress-energy tensor, 
and adopting a polytropic equation of state. The only 
nonzero components of the fluid four-velocity are u* and 
m''', and initial configurations can be constructed with 
either constant or non-constant specific angular momen- 
tum distributions, defined as j = hu^. We refer to [s^l 
for a full explanation of the construction of the initial 
model. For the simulations reported in this paper, we 
have considered a torus around a Schwarzschild BH (of 
mass Mbh = !)• The adiabatic index is F = 4/3 to 
mimic a degenerate relativistic electron gas, and the poly- 
tropic constant k is fixed to k = 0.0301262 such that the 
torus-to-BH mass ratio, Alt /Mbh, is roughly 0.1. 

The initial model is computed on a spherical grid. 
Therefore, in order to use it as initial data for our evolu- 
tion code, we transform into Cartesian coordinates and 
interpolate the values of the metric functions and hydro- 



dynamic variables onto a cell-centered Cartesian uniform 
grid, retaining data for the y = plane. 

The chosen torus has a constant distribution of spe- 
cific angular momentum such that its inner and outer 
edges on the equatorial plane are located at rjn — 7.1M 
and Tout = 14. OM, where M is ADM mass of system. 
Thus, such torus is initially covered with approximately 
140 points along the x-axis. The center of the torus is 
defined as the location at which the rest-mass density 
reaches its maximum and it is located at r,„ax = 9-8M. 
The dynamical timescale, which we choose as the orbital 
period at the center of the torus, is torb — 223M, that 
corresponds to about 1.1ms for the case that M — Mq. 

For these simulations we use an equally spaced uniform 
grid with a grid spacing Ax = Az = 0.05M and a grid 
with Nx xNz = 600 x 600 points to cover a computational 
domain, < a; < L and < z < L, with L = 30M. We 
use 4th-order finite differencing for the spacetime evo- 
lution, and Roe solver with MC reconstruction for the 
hydrodynamic evolution. The time integration is done 
with RK4. 



B. Dynamics in a fixed spacetime 

We have first investigated the equilibrium of the torus 
by performing numerical evolutions in a fixed spacetime, 
which nevertheless has the contribution coming from the 
self-gravity of the torus itself. This permits for one more 
test of the ability of the code to keep the stationarity of a 
fluid configuration initially in equilibrium. The accuracy 
check of the evolution is performed by comparing the 
stationarity and conservation of different local and global 
fluid quantities over a timescale which is several times the 
dynamical one. 

We show in Fig. [TO] the isocontours of the logarithm 
of the rest-mass density of the torus as computed at 
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FIG. 11: Time evolution of the total rest-mass, central rest- 
mass density and total angular momentum, each of them 
normalized to its initial value, for the evolution of a self- 
gravitating torus on a fixed spacetime. 



the initial time t = (left panel) and a,t t = llOOM 
(right panel), the later corresponding to 5 dynamical 
timescales, when the code was stopped after approxi- 
mately 3 X 10'' iterations, with no sign of the presence of 
numerical instabihties. These snapshots of the rest-mass 
density distribution clearly show the stationarity of the 
initial model for a sufficiently long period of time, and 
that the morphology of the torus remains unchanged for 
the duration of the simulation. In addition, they reflect 
that the treatment of the atmosphere in the vacuum re- 
gion does not perturb the equilibrium, and does not affect 
the dynamics during the evolution. 

Other quantities that provide a more precise measure 
of the stationarity of the torus are displayed in Fig. [TTl 
in which we plot the time evolution of three fluid quanti- 
ties normalized to their initial values. On the top panel, 
we plot the time evolution of the total rest-mass of the 
torus, normalized to its initial value. At the end of the 
simulation, the difference between the final total rest- 
mass with respect to the initial value is less than 10~^. 
This decrease of 0.1% of the total rest-mass is due to 
numerical error during the evolution, as the rest-mass 
is not conserved exactly. However, this error does not 
affect significantly the stability of the torus for the du- 
ration of our simulation. Similar results are obtained for 
the time evolution of the other two quantities, the cen- 
tral rest-mass density and the total angular momentum, 
which are displayed in the lower two panels of Fig. Ill[ 
respectively. The central rest-mass density, after a short 
initial transient phase, settles down to a stationary value 



FIG. 12: The upper panel shows the time evolution of the nor- 
malized central rest-mass density of a perturbed torus model, 
while the lower panel displays the power spectrum obtained 
after Fourier transforming the time evolution of the central 
rest-mass density. The fundamental mode and first overtone 
frequencies are denoted by / and oi. The mass of the BH is 
taken to be IMq. 



which only differs after t = lOOOM by 1% from the ini- 
tial one. This provides a strong evidence of the ability 
of the code to keep the torus in equilibrium for evolu- 
tions longer than the characteristic dynamical timescales 
of these objects. 

Notwithstanding the use of Cowling approximation, we 
can nevertheless extract some information about the os- 
cillating behavior of these type of self-gravitating tori. 
We note that in a series of recent papers, [zl, O, I?!, 
[tgI . \7T\ , it was shown that upon the introduction of per- 
turbations, non self-gravitating relativistic tori in equi- 
librium manifest a long-term oscillatory behavior lasting 
for tens of orbital periods. An important feature of these 
axisymmetric p-mode oscillations of accretion tori is that 
the lowest-order eigenfrequencies appear in the harmonic 
sequence 2 : 3. Overall, it was found that the 2 : 3 har- 
monic sequence was present with a variance of ^ 10% 
for tori with a constant distribution of specific angular 
momentum and with a variance of ^ 20% for tori with 
a power-law distribution of specific angular momentum. 
The departure from the 2 : 3 harmonic sequence depends 
on a number of different elements that contribute to small 
deviations, such as the vertical size of the tori, the BH 
spin, the distribution of specific angular momentum, the 
EOS considered, and the presence of a small but nonzero 
mass-loss, which can all influence this departure. 

In order to trigger the phase of small oscillations of the 
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FIG. 13: Isocontours of the logarithm of the rest-mass density 
of the torus during the evolution on a dynamical spacetime 
at t = 600M. (Compare with the left panel of Fig. ITUI ) 
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self-gravitating torus in our simulations, we perturb the 
equilibrium solution by adding a small perturbation in 
the a;-component of the velocity (we recall that in equi- 
librium the X and z-components of the velocity are zero). 
The oscillations reflect in the time evolution of the dif- 
ferent fluid quantities, for instance, the central rest-mass 
density, which we Fourier-transform to obtain the result- 
ing power spectra. This shows distinctive peaks at the 
frequencies that can be identified with the quasi-normal 
modes of oscillation of the disk. In Fig. [1^ we present, 
in the upper panel, the time evolution of the normalized 
central rest-mass density of the perturbed model. Due 
to the initial perturbation, the torus shows a persistent 
oscillating phase around its equilibrium position. We fol- 
low the evolution for about 13 dynamical timescales, and 
compute the power spectrum obtained from the normal- 
ized central rest-mass density, which is shown in the lower 
panel. A rapid look at this figure reveals that the torus 
power spectrum shows features which are very similar to 
those mentioned above for non self-gravitating tori (see 
e.g. Fig. 4 of [zE|)- Namely, it can be clearly identi- 
fied in the spectrum a fundamental mode /, and a first 
overtone oi . Interestingly, the ratio of the fundamental 
mode and the first overtone also show approximately the 
2 : 3 harmonic relation with oi/f ~ 1.4. In addition, 
the power spectrum reveals a series of higher frequency 
peaks, which roughly coincide with linear combinations 
of the / mode and first overtones. However, simulations 
lasting for longer timescales are needed to unambiguously 
identify these peaks. 

C. Evolutions on a dynamical spacetime 

Next, we evolve the same torus and BH initial data in 
a fully dynamic way, solving the Einstein equations cou- 
pled to the general relativistic hydrodynamics equations. 



FIG. 14: Time evolution of the total rest-mass density, cen- 
tral rest-mass density and total angular momentum, each of 
them normalized to its initial value, for the evolution on a 
dynamical spacetime. 



For such simulations we use the gauge conditions given 
in Eqs. (|2.22p - (|2.24p together with the x- method, which 
allow for long-term and stable evolutions of the puncture. 

In our simulations of the system formed by a BH plus 
a self-gravitating torus, the spacetime evolution is highly 
dynamical until t ~ 30M due to the initial adjustment 
of the gauge, which produces a small pulse in the metric 
functions that propagates outwards. Outgoing radiative 
boundary conditions for the spacetime variables at the 
outer boundaries, permit this initial pulse to leave the 
grid. Despite this initial transient phase, the torus re- 
mains in equilibrium during the evolution, and its mor- 
phology is kept very close to the initial profile even sev- 
eral hundred M beyond the first orbit when we stopped 
our simulation. In Fig. I13[ we plot the isocontours of 
the logarithm of the rest-mass density of the torus at 
time t = 600M, which is close to 3 dynamical timescales. 
This figure clearly demonstrates that the configuration 
remains in equilibrium. Other quantities also exhibit this 
behavior. In particular, in Fig. [TH we plot the time evo- 
lution of the total rest-mass of the torus (top panel) , the 
central rest- mass density (middle panel), and the total 
angular momentum (lower panel) , each normalized to its 
initial value, for the duration of the simulation. At the 
end of the simulation, the difference between the final to- 
tal rest mass with respect to the initial value is less than 
0.1%. Although it is due to accumulated numerical error 
during the evolution, this violation of the rest mass does 
not affect significantly the stability of the torus for the 
duration of our simulation. Similar results are obtained 
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for the time evolution of the other two quantities, the 
central rest-mass density and the total angular momen- 
tum. We note that, truncation errors at the resolution 
employed seem to be enough to trigger small oscillations 
of the torus around its equilibrium as shown in the evo- 
lution of the central rest-mass density. 

It is worth to mention, that in order to maintain the 
torus in equilibrium for several hundred M, we find cru- 
cial not to replace the by —djj'^^ on the right hand side 
of the evolution Eq. (|2.9p whenever is not differentiated. 
Replacing F* by —djj^^ excites larger amplitude oscilla- 
tions of the torus, which can be more than 10% larger 
than in the case where F' is not replaced by its defini- 
tion. We note also that, as observed for single puncture 
simulations in vacuum, numerical errors in the spacetimc 
evolution, especially near the puncture, are smaller with 
the x-method than with the (/(-method. We find that 
an additional ingredient that helps to maintain the torus 
equilibrium for the duration of the simulations is the use 
of the x-method for the evolution of the conformal fac- 
tor (we refer to Fig. [3] and Fig. d] for a comparison of the 
X-method and 0- method). We find numerical indications 
that, at least in conjunction with the Cartoon method, 
these two procedures help to decrease the error in the 
violation of the Hamiltonian constraint. 



VII. CONCLUSION 

In this paper we have presented a new two-dimensional 
numerical code designed to solve the full Einstein equa- 
tions coupled to the general relativistic hydrodynamics 
equations. The code is mainly intended for studies of 
self- gravitating accretion disks around BHs, although it 
is also suitable for regular spacetimes. 

Concerning technical aspects, the Einstein equations 
are formulated and solved in the code using a reformu- 
lation of the standard 3-1-1 (ADM) system, the so-called 
BSSN approach, in conjunction with the Cartoon method 
to impose the axisymmetry condition under Cartesian co- 
ordinates, and the puncture/moving puncture approach 
to carry out BH evolutions. We note that a key and 
novel feature of the code is that combines a fourth-order 
finite difference approximation for derivative terms in the 
spacetime evolution with the Cartoon method. Corre- 
spondingly, the general relativistic hydrodynamics equa- 
tions are written in fiux-conservative form and solved 
with high-resolution, shock-capturing schemes. 

We have performed and discussed a number of tests to 
assess the accuracy and expected convergence of the code, 
namely (single) BH evolutions, shock tubes, and evolu- 
tions of both spherical and rotating relativistic stars. We 
have also presented a simulation of the gravitational col- 
lapse to a BH of a marginally stable spherical star. We 
have shown that the code is able to handle the forma- 
tion of a BH, and to follow the BH evolution until we 
the simulation is stopped after several hundred M. We 
remark that we did not add any numerical dissipation to 



the evolution equations for the spacetime variables and 
gauge quantities to perform this simulation; and exclu- 
sively relied on the gauge choice to follow the BH forma- 
tion and its long-term evolution. Overall, the code has 
passed those tests with remarkable accuracy. 

In addition, paving the way for specific applications of 
the code, we have also presented results from fully gen- 
eral relativistic numerical simulations of a system formed 
by a BH surrounded by a self-gravitating torus in equi- 
librium. First, simulations using a fixed spacetime evo- 
lution, have shown that the torus remains in equilibrium 
around its initial configuration for more than 5 dynamical 
timescales, when the simulation was stopped. Further- 
more, after adding a small perturbation on the torus in 
equilibrium, we have followed the evolution of the torus 
through a phase of small oscillations for more than 13 
dynamical timescales. The computation of the frequen- 
cies of the fundamental mode of oscillation and the first 
overtone revealed that these two frequencies appear close 
to the 2:3 harmonic relation observed for the case of non 
self-gravitating relativistic tori. 

Our simulation of the same model in a dynamical 
spacetime, showed that the code is able to keep the torus 
in equilibrium for the several hundred M that lasted the 
simulation. We have found numerically, that using the x~ 
method and not to replace the F* by —dj^'^^ on the right 
hand side of the evolution for the F* are two important 
issues that reduce the numerical error in the spacetime 
evolution, and therefore, help to maintain the torus in 
equilibrium. 

Overall, the simulations performed indicate that the 
code is able to perform such simulations accurately and 
for the sufficient duration needed to produce scientific re- 
sults for one of the specific applications for which it has 
been designed, the investigation of the runaway insta- 
bility of thick accretion disks around BHs. We point out 
that the last simulation presented in the paper is the first 
2D axisymmetric simulation of the system formed by self- 
gravitating torus around a BH in a dynamical spacetime 
that has been ever carried out. In a follow up paper, we 
aim to present results from an expanded range of initial 
models, and to investigate in detail the dynamics of such 
systems, focusing on the development on the runaway in- 
stability, and the evolution of the BH due to the transfer 
of mass and angular momentum from the torus. 
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